import numpy as np

from pynumpak.type import func_t, Interval


def romberg(f: func_t, ab: Interval, tolerance: float = 1e-7):
    """romberg 法求积分"""
    a, b = ab
    interval = b - a
    last_row = [(f(a) + f(b)) * interval / 2]

    for i in range(1, 100):
        points = np.linspace(a + interval * 2 ** (-i), b - interval * 2 ** (-i), num=2 ** (i - 1))
        row = [last_row[0] / 2 + np.sum(f(points)) * interval * 2 ** (-i)]
        for j in range(i):
            result = (4 ** (j + 1) * row[j] - last_row[j]) / (4 ** (j + 1) - 1)
            row.append(result)

        if abs(last_row[-1] - row[-1]) < tolerance:
            break
        last_row = row
    return result
